Quantum State Tomography: 'the best' is the enemy of 'good enough'. 
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In this paper, we examine a variety of strategies for numerical quantum-state estimation from 
data of the sort commonly measured in experiments involving quantum state tomography. We find 
that, in some important circumstances, an elaborate and time-consuming numerical optimization to 
obtain 'the best' density matrix corresponding to a given data set is not necessary, and that cruder, 
faster numerical techniques may well be 'good enough'. 
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I. INTRODUCTION 

The goal of quantum state tomography [11 [H |3] is to 
estimate, from a series of projective measurements per- 
formed on identically prepared quantum systems, the 
density matrix of the underlying ensemble of which these 
quantum systems are realizations. This process is nec- 
essarily non-deterministic in nature, relying on the fre- 
quency of experimental outcomes to estimate probabili- 
ties - a process that converges to the actual probabilities 
only in the infinite limit. Thus the reconstruction of the 
quantum state cannot be exact in any realistic experi- 
ment. Furthermore, these measurements can only yield 
estimates of the on-diagonal elements of the density ma- 
trix, but not directly any data about the off-diagonal 
elements. It is necessary to perform various unitary op- 
erations on the system (or, equivalently to perform pro- 
jective measurements in a variety of bases) in order to 
obtain such information about the complete state. In- 
deed, for a system with a discrete spectrum of n-levels, 
the density matrix is specified by — 1 independent 
real parameters, and each parameter will require a sep- 
arate measurement. Even after the required measure- 
ments have been performed, the experimenter faces the 
problem of estimating the density matrix from incom- 
plete and noisy data. The problem is aggravated by the 
constraints that quantum physics places on the density 
matrix: It must be a non-negative, unit-trace Hermitian 
matrix. Today, the approach that is usually taken is to 
determine computationally what is the 'best' such pos- 
itive, unit-trace Hermitian matrix which corresponds to 
a particular data set, and what confidence can we place 
on such an estimate. The most complicated such to- 
mographic measurement performed to date on an 8 
qubit (256 state) system, realized in a trapped ion exper- 
iment, was limited not by the experimental capabilities of 
the system, but rather by the complexity of the numeri- 
cal state recovery problem [B]. This computational com- 
plexity, while underscoring the awesome computational 
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potential inherent in quantum information, nevertheless 
presents an experimenter, intent on exploring larger and 
larger Hilbert spaces, with considerable tribulation when 
characterising the performance of his or her apparatus. 

In this paper we examine the problem from an entirely 
computational perspective. Specifically, we address the 
concern that maybe we are being too fastidious in ap- 
proaching the state reconstruction problem. One can 
obtain a positive, unit-trace Hermitian matrix from to- 
mographic data in a variety of ways. First, and most 
simply, one could generate a linear reconstruction of the 
noisy data (which tends to give a non-positive matrix), 
and ensure positivity by setting the negative eigenvalues 
to zero, then re-normalizing to ensure a unit-trace. This 
we call the "Quick and Dirty" (QD) approach. A second 
strategy is to assume the state must be nearly pure - after 
all, quantum technologies are usually in the business of 
trying to create pure states - and to simplify the compu- 
tation by finding the pure state most compatible with the 
data. We call this the "Forced Purity" (FP) approach. 
A third approach is full optimization, i.e. the application 
of some constrained optimization routine, with a specific 
metric to define the 'distance' between our data set and 
a positive density matrix, and search parameter space 
until the absolute 'best' (i.e. global minimum) density 
matrix is obtained. Our goal is specifically to address 
the question: When is the rigorous optimization required, 
and when will some short-cut technique be good enough? 
This is a question that can only be addressed by simu- 
lation: Since we need to know a priori the underlying 
density matrix of the ensemble to compare the recovered 
estimates. Starting with assumed density matrix, we em- 
ploy a pseudo-random number generator to create some 
'pseudo-experimental data' with appropriate probability 
distribution. The various approaches to density matrix 
recovery are applied to it, and the result is compared 
with the initial density matrix to assess the accuracy 
of the recovery techniques. Our analysis concerns solely 
multiple correlated two-level systems, e.g. the qubits of 
a small scale quantum computer; however, many of the 
techniques and results we present are readily adaptable 
to more general systems. 

The paper is organized as follows: In Sec. [H] we discuss 
the generic tomography problem for a single qubit, which 
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is generalized to the n-qubit case in Sec. Ill describing 
specifically a number of memory management techniques 
required for scalability of the code, and our novel ap- 
proach to the optimization routine (using gradient-based 
algorithms and employing the matrix differential calcu- 
lus). The code itself is described in detail in Sec. IV and 
our results in Sec. IVl 



II. ONE QUBIT 

In this section, we will review the basic concept of 
quantum state tomography by considering the estima- 
tion of a state of a single two-level system, or qubit. 



A. Parametrizing the Density Matrix 

The density operator describing the state of a sys- 
tem [7] is a Hermitian, non-negative definite operator of 
unit-trace. The set of Pauli matrices [S] {ctq, (Ti, i5'2, ^g} 
form, for a two dimensional space, a complete orthonor- 
mal set of matrices, so that p can be expanded as a linear 
combination of a,, as 



(1) 



where 



= Tr{a,p}/2. (2) 

Since Tr{p} — X^r^ — 1/2; further, since — p. the r,y 
are all real parameters. 

The may be determined experimentally as follows: 
Suppose we perform a measurement, specified by the pro- 
jector Ho, on the system, the probability of obtaining a 
positive outcome is TrjpIIo}. Repeating this measure- 
ment M times on identically prepared systems, the ex- 
pected number of times we obtain this outcome will be 



no = AfTr{ptl^} 
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(3) 



If one repeated this procedure of multiple measure- 
ments for a set of four different measurement operators, 
{III/}, (i' — 0, 1, 2, 3) one obtains a set of Hnear equations 



where 



(4) 



(5) 



By choosing the measurement operators, {ft^}, judi- 
ciously, one can ensure that B^,^^ is non-singular, and 



hence that the desired parameters can be obtained 
from the observed quantities n^, viz., 



(6) 



M=0 



Substituting r^, into Eq. Q, we obtain the density ma- 
trix, as a function of measurement outcomes, provided 
the measurements have no noise or errors in them. 

Following the precedent of Ref. ^ we use the standard 
Stokes measurement basis for our numerical experiments. 
These measurement operators are given by: 



no = i(|o)(o| + |i)(i|), ni = 

Il2 = \D){D\, U3^\R){R\, 



(7) 



where |0) and |1) represent the two states of our qubits, 
and 



\R) 
\D) 



1 

1 



(|o)-»-|i)), 
(|o)-|i)). 



(8) 
(9) 



A natural metric to compare the recovered density ma- 
trix Pmeas with the actual density matrix ptruc is the fi- 
delity [19], defined as: 



F{p 

mcas 7 Ptruc mcasPtrucV Pmcas 

(10) 

However, when we invert the measurement data linearly, 
our recovered "density matrix" piincar is not non-negative 
definite and hence we have the specific problem that fi- 
delity turns out to be complex (not to mention the more 
general problem that piinear cannot be interpreted as a 
density matrix of a physical state). We have to correct 
the matrix obtained by linear reconstruction to obtain a 
proper density matrix. 



B. "Quick and Dirty" Reconstruction 

As a simple initial approach to this problem, we can 
decompose puncai into its spectral representation, i.e. 



Piii 



(11) 



where D is the diagonal matrix of eigenvalues (which 
are real, but not necessarily positive) and [/ is a unitary 
matrix. We then set all negative eigenvalues in D to zero, 
call this matrix D' and obtain: 



PQD 



Tr{D'}' 



This provides a rough initial estimate of the state; one 
of the goals of our analysis in this paper is to assess how 
good an estimate it is. 
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C. "Forced Purity" 

An alternative approach to the problem of obtaining a 
non-negative definite density matrix from measured data 
is to assume that the state is pure. Recall that for a pure 
state \^) the density matrix can be described by a single 
ket as Ppurc — Such a density matrix for n qubits 

has eigenvalue with degeneracy 2" — 1 and eigenvalue 
1 with degeneracy 1. 

Because ppurc is also Hermitian, it can be written in 
its spectral decomposition as 



Ppu 



where impure is the diagonal matrix with a single element 
equal to 1, and all other elements being zero; is a 
unitary matrix. 

During linear inversion of a pure state, the eigenvalues 
of piincar may be negative, but sufficiently close to eigen- 
values of Ppure- The idea of forcing purity on such a state 
is to obtain 



Pfp 



1^ linear^ ''linear 

Tr{z5"} 



where D" is the diagonal matrix obtained from D by 
setting the largest eigenvalue equal to 1, and all others 
equal to 0. 



D. Maximum Likelihood 

Any Hermitian 2x2 non- negative unit-trace matrix can 
be uniquely parametrized using the Cholesky decompo- 
sition as: 



Pidcal {tl, t2, ^3, ^4) 



Tr{T^T} ' 



where 



T(ti, t2, ^3, t4 



tl 



(12) 



(13) 



Thus a 'physical' density matrix can be specified by the 
four parameters t = {^i, t2^ ts, t^}. The ideal of the maxi- 
mum likelihood method is to perform a search of the t pa- 
rameter space until we find a /Oidcai (i) which is most likely 
to have generated the observed data {no, ni, 7x2, 713}. To 
assess this likelihood, suppose that each datum is 
a statistically independent, Poisson-distributed random 
variable with expectation value h^. Further, if is a 
large number, the Poisson distribution is well approxi- 
mated by the Gaussian distribution, i.e. 



P(no,ni,n2,n3) 



1 



n 



1^=0 



(n^ ~ n ^) 
2n, 
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(14) 



each datum is garnered from J\f repetitions of 
a measurement carried out on a system in state 
Pideai {t) , it is reasonable to make the identification 
n„{ti,t2,t3,t4) = J\f{ip„\pidca.i(ti,t2,t3,U)\4'iy), and the 
likelihood of a given parameter vector t generating the 
data {no, n-i, ^2, ^,3} can be obtained by substituting this 
identity into Eq. ( 14 1 . We are then in a position to de- 



termine the parameter vector for which this probabil- 
ity is maximized, and hence the most likely density ma- 
trix. Instead of maximizing Eq. ( 14 1, it is equivalent, and 



mathematically more convenient to minimize the follow- 
ing function: 



1 3 ^n^-7Vrr{n 

2 ^ 



i/Pideal 



(15) 



In order to optimize this function efficiently, we need 
to compute its gradient. This is not an easy feat, as the 
closed analytic form does not simplify well, and finite- 
differencing is too inefficient. The situation becomes ex- 
ponentially worse as we increase the number of qubits. 



III. GENERALIZATION TO N-QUBITS 

In the previous section, we outlined the possible rou- 
tines for performing tomography of a single qubit. We 
now extend these routines to a higher number of qubits 
and see how the "Quick and Dirty" and "Forced Purity" 
methods compare to the elaborate and time-consuming 
Maximum Likelihood Estimation (MLE) routine. 

At first, the problem looks very simple - any state of 
each qubit is completely characterized by only 4 measure- 
ments. Hence, numerically the MLE procedure is rather 
easy to implement - we just need to optimize a func- 
tion of 4 variables, which is achieved by the simplex or 
Powell optimization algorithm in a fairly short amount 
of time [9 , without computing the gradient. However, 
two qubits, when correlated, are not characterized by 
8 measurements, but by 4 • 4 = 16 measurements, be- 
cause we are looking at a system of 2 qubits. If n is the 
number of qubits, then we would need to obtain 4" mea- 
surement outcomes in some fixed 4" dimensional basis. 
Due to wave function collapse, we can only perform one 
projection measurement at a time (an outcome is an aver- 
age over multiple identical projection measurements) and 
for each projection measurement on one qubit, we have 
to cycle through all possible combinations of projection 
measurements for the other qubits. 

Let us introduce the following set of operators which 
generalize the Pauli matrices for n qubit systems: 



1 



/2" 



(16) 



where < /i^ < 3 for all 1 < ^ < n arc the digits of 
the index p in base-4. For example, if /i = 33 for a 4- 
qubit system, r33 = (To (E) o'2 (E) (To (8) (Ti , since 33 is equal 



4 



to 0201 in base-4. For convenience, we have included 
a normalization constant, so that Tr{TfjT,^} = 5^^^ (in 
keeping with the convention used in Ref. |3]). Similarly, 
we write the projection operators for our measurement 
states as 



= (g) Hj,, 



(17) 



The Cholesky decomposition of p remains the same, ex- 
cept that T{t} is a 2" X 2" matrix specified by 4" param- 
eters t — {ti,t2, . ■ ■ ^4"}, i.e. 



h 

^2"+l + i • ^2"+2 ^2 
^4"-l + i - tin. 



<2"+i-4 + * • ^2"+i-3 ^2" 
(18 



A. Computational Constraints and 
Memory-efficient Linear Reconstruction 

In order to perform computational numerical tomog- 
raphy in practice, we need to take the following into con- 
sideration: 



1. Computational efficiency - what is the upper bound 
on the number of ffoating point operations of a cer- 
tain tomography algorithm. 



2. Amount of memory available - what is the up- 
per bound on the size in computer memory of the 
largest data structure used by the tomography al- 
gorithm. 



Kronecker tensor products increase the size of resultant 
matrices exponentially. The goal is to obtain the den- 
sity matrix which has 2" x 2" elements. So we cannot 
have any other data structure in memory which would be 
larger, otherwise the problem of increasing the number 
of qubits becomes constrained by that particular data 
structure. 

For example, consider the approach described in 
Ref. |4j, in which a 4" x 4" complex matrix -B^.jy (the n- 
qubit generalization of the matrix defined by Eq. ^) was 
stored in memory. The table below outlines how much 
memory is needed to store a 4" x 4" complex ffoating 
point matrix, using 32 bits to store the real or imaginary 
part: 



Qubits 


Bytes 


GigaBytes 


1 


128 


1.28 X 10"^ 


2 


2.05 x 10^ 


2.05 X 10"® 


3 


3.28 x lO"* 


3.28 X 10"^ 


4 


5.24 X 10^ 


5.24 X 10"" 


5 


8.39 X lO'' 


0.01 


6 


1.34 X 10^ 


0.13 


7 


2.15 X 10^ 


2.15 


8 


3.44 X 10^" 


34.4 


9 


5.50 X 10" 


580 


10 


8.80 X 10^^ 


8.80 X 10^ 


11 


1.41 X 10" 


1.41 X 10^ 



Table I: Amount of memory required to store a 4" x 4" com- 
plex floating point matrix using 32 bits to store the real or 
imaginary part. 

It must also be noted that any type of storage media 
has to be able to perform read and write operations quite 
fast because this data structure would be accessed quite 
frequently. This is simply not the case for most conven- 
tional hard drives: Using standard personal computers 
of the type typically integrated into quantum optics lab- 
oratories, one is in practice limited to about 7 qubits, 
without resorting to more powerful computer hardware. 
However, a data structure of maximum size of 2" x 2" 
would allow to go as high as 15-16 qubits, at which point 
the density matrix itself would become a storage prob- 
lem. Thus our goal is to avoid storing matrix -B^.^ into 
memory. Instead, we can obtain its inverse element by 
clement. This can be achieved as follows: The matrix 
B^i^i, for an rt-qubit system is defined by the equation 

= Tr{n,,a^jrr{n,,a^J...rr{n,„a^J.(19) 

Defining the 4x4 matrix f3v^,tL^ ~ Tr\tl^^a^^} for all 
1 5: C 5: '^i which can be easily inverted (provided a 
suitable set of measurements {11/^}, (/i = 0,1,2,3) has 
been chosen) , we find 



(20) 



where as before, (z^i, f2, ...z^n) are the base-4 digits of the 
index u (and similarly for /i). 

This allows us to calculate the initial linear recon- 
struction of the density matrix from the observation data 
(no,ni, ....Ui^^i), viz: 



where 



Piii 



(AA) 



4"-l 

E 



(21) 



(22) 



5 



in a computationally efficient manner. Now the only size 
constraint on linear reconstruction is the density matrix 
itself. Of course, storing the projection measurement ma- 
trices III, is also problematic - a quick solution is to gen- 
erate these matrices when they become needed - one can 
store certain tensor combinations which make up 11^ into 
memory and only tensor on additional combinations to 
obtain the desired llj/. 



B. Maximum Likelihood 

Extending the Maximum Likelihood Function (MLF) 
from Eq. (151 to n qubits we obtain: 



4"-l 



[AArr{n,pidoai(i)} 



(23) 



where to simplify calculations we assumed that we can 
approximate variance by the measurement outcome av- 
erage in the denominator. Minimizing this function be- 
comes a severe computational problem. Most gradient- 
free optimization routines are rather slow and only work 
well for a low number of dimensions, whereas here we 
have a number of dimensions which grows exponentially 
with the number of qubits. Wc have to use a numeri- 
cal routine which is more efficient; this usually involves 
calculating the gradient and/or the Jacobian. The finite- 
differencing approach is too slow for computing the gra- 
dient (a fact which we verified computationally) because 
evaluating the MLE function is exponentially inefficient. 
Hence, we require an analytic closed form for the gradient 
and/or the Jacobian matrix. 

It should also be noted that if the region of opti- 
mization is convex, we are looking at non-linear convex 
optimization problem, for which a number of algorith- 
mic approaches should work. We decided to take the 
simplest approach possible: Optimize the MLE function 
with built-in constraints using an algorithm which works 
on both convex and non-convex sets. We reduce the com- 
putation time by deriving an analytic form for the gradi- 
ent. An alternative approach is to derive a different MLE 
function with an external set of constraints, and launch 
another convex optimization algorithm similar to linear 
programming [15i il6j . 



1. Initial Algorithmic Attempts 

The following algorithms were considered to optimize 
C [9l [10], mostly because they are available in libraries 
such as GNU Scientific Library (GSL) [20]: 

1. simplex method; 

2. Powell's quadratically convergent method; 

3. Levenberg-Marquardt nonlinear least squares; 



4. conjugate-gradient method; 

5. BFGS algorithm. 

Routines |4] and |5] need to be able to perform gradient 
computation and line search in an efficient manner and 
have to converge to the desired minimum, even if started 
far away from it. 

Let us start with the line search routines - the following 
algorithms can be implemented for line search: 

1. successive parabolic interpolation; 

2. Newton's method; 

3. Golden Section Search (GSS). 

In order to pick one algorithm out of these three, we need 
to first know if the region of optimization is convex or not, 
and if so then how closely does it resemble a quadratic 
function. 



2. Matrix Calculus Derivation 

Regardless of the method chosen for the line search in 
Sec. |IIIB 1| wc still need an efficient way of computing 
the gradient. As established earlier, finite-differencing re- 
quires too many function evaluations and does not satisfy 
our computational efficiency constraints. 

The goal of this section is to find the gradient of C 
in closed form. This procedure can then be extended to 
finding second order partial derivatives for the Hessian 
matrix and differentiation with respect to a constant for 
line search routine. 

We begin with the gradient derivation. This reduces 
to finding 



9£ 



dC dT 



(24) 



Using matrix calculus, it suffices to find which would 
be a matrix of size 2" x 2" in our case. Certain elements 
of this matrix would represent the values of ^: 



dT " Tr{rt(i)r(i)}2 



E 



7Vrr{H^Pidcai(i)} - n^, 



Tr{r^(t)T(i)} 



arr{n,Tt(i)r(t)} 



dT 



dTr{T^i)T{t)} 



dT 



(25) 



Defining the following real quantities: 



A{?) = Tr{T^i)T{i)} and 
B,(i) = Tr{U,T^P)T{t)}, 



(26) 
(27) 
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we find: 



Tr{fl^Pidcai{i)} 



A 



(28) 



Further we will denote the matrix derivatives of these 
quantities with respect to the Cholesky matrix T as fol- 
lows: 



A'{T) = 



dTr{fl,T^^)T{t)} 
dT 

dTr{T^i)T{^} 
df ■ 



(29) 



(30) 



Because matrix calculus is only well-defined for real- 
valued matrices, let us write 



T = T{t) =X + i-Y, 



(31) 



Then, using the matrix calculus theorems in Sec. \K\ we 
find 

A' = 2X + i-2Y, (32) 
Bl = 2XK^-2YA^ + i-(2XA^ + 2YK^). (33) 



Denoting 



NBt, - An^ 
An,, 



ABl - B,A' 



A^ 



(34) 



we find the matrix derivative of C can be written in the 
compact form 



dC 

df 



£'(*)- CD,, 



(35) 



where is a scalar and is a 2" x 2" matrix. In 
fact, the upper-diagonal of C'{i) and imaginary part of 
the diagonal are of no use to us - values of the gradient 
are seeded in the original locations of t,, so C'{i) has to 
be disassembled into real and imaginary parts and then 
the gradient vector has to be filled from the resulting 
matrices. 



IV. DESCRIPTION OF CODE 

The goal is to scale tomography routines up to a higher 
number of qubits on a standard single-processor worksta- 
tion by refining the tomography algorithms to remove the 
numerical complexity bottleneck from experimental post- 
processing. In this section, we describe how the codes 
were implemented. 



1. Linear reconstruction - provides the linear recon- 
struction of the data by inverting the measurements 
into a matrix puncar^ outlined in Sec. [Ill A[ which 
has all characteristics of a density matrix, except 
positive semi-definiteness. 

2. "Quick and Dirty" - quickly fixes piinoar into /5qd 
by setting all negative eigenvalues to zero and re- 
normalizing. 

3. "Forced Purity" - for pure states, eigenvalues are 
known. This routine forces eigenvalues of piincar or 
Pqd (does not matter which one) into those of a 
pure state, also ensuring unit-trace condition. 

4. MLE - we use the elements of the "Quick and 
Dirty" density matrix as a starting point for our 
optimization routine. We then launch the BFGS2 
algorithm supplied with GSL. 

Our progress while developing these routines is as follows: 

1 . Started with our own simplex method code in Mat- 
lab, which only optimized 4 qubits - gradient-based 
algorithm was needed. 

2. Wrote the conjugate-gradient routine in Matlab 
using GSS line search routine and using finite- 
difference gradient, which allowed for 5-6 qubit to- 
mography. 

3. Applied matrix differential calculus to the gradi- 
ent and obtained a closed form expression, which 
severely improved Matlab routines for up to 7 
qubits. 

4. Experimented with Newton's method and succes- 
sive parabolic interpolation line searches which did 
not work in the end. This led us to suspect that 
the region of optimization is not convex. 

5. Re- wrote everything in C using GSL and employed 
GSL's BFGS2 algorithm and its collection of line 
searches - this pushed our routines to 9 qubits 
(MLE hmits to 9 qubits, but not "Forced Purity"). 

All code is currently implemented in C using GSL, with 
prototype routines also available in Matlab. 



B. Creating Pseudo-Experimental Data 

Generally, if one wants to simulate a physical state 
characterized by /Sphysicai with 100e% experimental state 
error, then 



A. State Tomography Routine 



Pphysical — Pthcorctical ~t- ' ^random; 



Four routines provide tomography and run in the fol- 
lowing order: 



where Ptheoreticai IS & density matrix of some desired state 
and e, a real- valued constant, simulates experimental 
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"state error" - the physical state always differs from the 
intended state by some small amount; a random density 
matrix is created as follows: 

i? = 2 • rand(2'") - I + i ■ (2 ■ rand(2") - 1), (36) 



/^random " 



Tr{R^R}' 



(37) 



where rand function creates a 2" x 2" matrix of pseudo- 
random values, sampled from Uniform(0,l) distribution 

For instance, the following results in a noisy GHZ state: 

PGHZ - (1 - e) • ^|100...01)(100...01| + e • Prandom- 

The simulation routine creates a physical density matrix, 
simulates experimental measurement outcomes and then 
attempts to reconstruct this density matrix. Knowing 
what the reconstructed density matrix should be, enables 
us to compare how well each reconstruction routine works 
for a certain number of qubits. 

The expected number of positive outcomes is obtained 
usmg Eq. ^, viz: 



' n^/^ ^physical} ; 



where A/" is a constant which is equivalent to the number 
of times repeated projective measurements were taken 
|22j . We then add experimental noise to the measure- 
ments [53] using: 

= Poisson{n^), 

where Poisson{X) generates a random number from a 
Poisson distribution with mean A using a probability in- 
tegral transformation. 



for n qubits. 

The tangle (i.e. the square of the concurrence [T7]) is 
defined for 2 qubits, as 

r = [max{\i - Ai - A2 - A3, 0}]^, 

where A's are the square roots of the eigenvalues of the 
matrix \/ p(py ® ay)p*{uy ® (jy)^, which is guaranteed 
to be Hermitian [17 and cTy ® Uy is the spin-flip matrix 
and p* is the complex conjugate of density matrix p. For 
larger numbers of qubits, it can be used as a lower bound 
on the degree of entanglement [T5]. 



A. Linear Entropy vs Tangle Plane 

We observed that for 2 qubits, certain states produce 
fidelities of over 90% using "Quick and Dirty" routine 
and consistently high fidelities using MLE. We generated 
2 X 10^ pseudo-random density matrices, that filled the 
entire entropy-tangle plane. Random density matrices 
had to be biased in order to evenly cover the entire plane. 
For example, to fill the plane below the Werner state line, 
we used: 



Pt = 



1-52 sVT^ 



SVl^ 



(38) 



which biased the tangle in 



Ptrial = • Prandom + (1 ~ C^) ' Pr, 



where 



V. RESULTS 

In this section, we discuss the conclusions to be drawn 
from the numerical trials described in the previous sec- 
tions. In particular we address the question posed in 
the title of this paper: Do we always need an expensive 
MLE routine to perform tomography or would "Quick 
and Dirty" or "Forced Purity" methods suffice? 

We compare the "Quick and Dirty" and "Forced Pu- 
rity" routines to the "MLE" routine for states with wide 
variations of entropy and entanglement. We also show 
how well these routines scale in runtime and how ex- 
perimental errors affect the reconstructed states as the 
number of qubits increases. 

The linear entropy, which specifies the degree of purity 
of the state, is defined as 

'S'li„ear(p) = ^„ _ [1 - Tr{p}] 



< e < 1, 



< (5 < 1 



Varying S changes tangle and varying e changes en- 
tropy. We cycled through 100 x 100 different combina- 
tions of S and e and for each setting performed 100 trials, 
sampling Prandom from Uniform(-l,l) probability distri- 
bution for each trial. We can also move along the Max- 
imally Entangled Mixed State (MEMS) line by varying 
7, 



Pmems 



5(7) 7/2 

1 - 23(7) 



7/2 3(7) 



where 



5(7) = 



7/2, 7 > 2/3 
1/3, 7 < 2/3 
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and 



Ptrial 



^ ■ Prandom + (1 ^ E^) ' PmEMS- 



We further sampled 1000 x 1000 different settings of e 
and 7 to fill the area around MEMS line: In this case 
increasing e increases the distance from the MEMS line. 




Figure 1: On the Suncar and r plane you can see that the 
2 X 10^ generated states covered the entire plane; above each 
point is the corresponding fidelity of the recovered state using 
maximum likelihood. Notice that most fidelity values lie be- 
tween 90%-99%. The fidelity values create a thin plane, which 
suggests that the standard deviation is low for various states. 
However, some high-entropy states cannot be recovered well 
even with the expensive MLE procedure. 



regardless of the value of t, share one thing in common: 
eigenvalues. More precisely, for n qubits, eigenvalue 
occurs with degeneracy 2" — 1 and eigenvalue 1 occurs 
with degeneracy 1. So, setting negative eigenvalues to 
zero adjusts the eigenvalues closer to the eigenvalues of 
a pure state. If we know that the state is pure ahead 
of time, we can just reset the eigenvalues to the known 
values after the linear inversion procedure and obtain the 
density matrix - this is further explored in Sec. V C 




Figure 3; Again, same 2 x 10 states, but this time with 
"Forced Purity" performed. Notice that for 2 qubits "Forced 
Purity" appears to be worse than "Quick and Dirty" even for 
pure states - this is not the case when the number of qubits 
increases, as we explore in the next section. 




B. Performance for N qubits 

In order to extend this assessment to larger numbers 
of qubits, while still varying the amount of entangle- 
ment and disorder, we considered a generalized version 
of the Werner state for n qubits. Since this state slices 
through the entire plane presented in Sec. V A we can 



see how well tomography operates on states with various 
tangle and entropy values by varying the location along 
the Werner state line. An adjusted Werner state density 
matrix is given by: 

PWcrncr = | GHZ) (GHZ | • 6 + til . ^ 



Figure 2: These are the same 2 x 10^ states as in the previous 
figure and the projection on Suncar and r plane is identical. 
Notice how the "Quick and Dirty" routine also fails to reach 
high fidelity values as the states become more mixed. How- 
ever, for pure states Quick and Dirty is comparable to MLE 
in fidelity values. 

This suggests that for states with low entropies (pure 
states). Quick and Dirty routine should work in theory. 
This is not surprising, as from the spectral decomposi- 
tion, we can see that all states with S'lincar = property. 



where 



|GHZ)--^{|00...0) + |11...1)} 



and / is the 2" x 2" identity operator, representing a 
maximally mixed state. When e = 1 we obtain a state 
located at Sunear — and r = 1 and when e = we 
obtain S'linear = 1 ^ud r = 0. We then vary e from 
to 1 in 101 increments and for each value of e perform 
tomography 100 times, for a fixed number of qubits. 
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Werner State parameter, e 



Figure 4: MLE works well for 2 qubits, but for a higher num- 
ber of qubits there is a significant drop in fidelity. A 10% state 
error is quite large, but even at this error pure states are re- 
constructed better than mixed states as you can see around 
e = 1 (highly entangled pure state). 



C. "Forced Purity" Tomography for Pure States 

In order for "Forced Purity" to work, the measurement 
outcomes have to be sufficiently close to their true val- 
ues. To address the issue of 'how close', we simulated 
pure states, with 11 distinct tangle values evenly spaced 
between and 1 and then started with A/" = 10 - i.e. A/" 
repeated projective measurements for each measurement 
outcome. We then increased Af by one at each itera- 
tion and repeated "Forced Purity" tomography for some 
number of qubits. As soon as Af allowed "Forced Purity" 
to perform tomography at 90% fidelity, the routine was 
terminated and the value for Af recorded. 

One of the possible reasons why MLE tomography did 
not yield high fidelity values for a larger number of qubits 
is that it also required more accurate estimates of the 
measurement outcomes. Because MLE produced almost 
equal fidelities to "Forced Purity" for pure states, we 
would expect the same number of Af to work for MLE 
tomography. 



m 0.5- 




2 qubits 

3 qubils 

4 qubits 
5 qubils 

- - 6 qubits 



Werner State parameter, i 



0.9 



1 



Figure 5: "Quick and Dirty" routine is not performing very 
well even for pure states, although for a certain number of 
qubits it appears to work. On the contrary, it appears to be 
improving as the number of qubits is increasing. 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Werner State parameter, e 



Figure 6: For "Forced Purity" we can confidently say that the 
routine is improving significantly, as the number of qubits 
is increasing. And for pure entangled states this routine is 
almost as good as MLE, but runs at a fraction of the time! 



D. Runtime Analysis 



Section V C suggests that pure states do not require 
expensive MLE techniques for tomography. Nonetheless, 
it is interesting to see how MLE scales in runtime com- 
pared to Quick and Dirty and "Forced Purity" routines. 
Here we present runtimes and fidelity estimates for a pure 
state with r « 0.5 and a slightly mixed state with the 
same tangle value which lies on the Werner state line. 
We also show that even the expensive MLE routine de- 
creases in fidelity as we increase the number of qubits. 
For this analysis we assume that an experiment can be 
performed a sufficiently large number of times, 10^ to be 
exact. 

In conclusion, "Forced Purity" results in lower Fidelity 
values for 2 qubits than "Quick and Dirty" , but then 
increases in Fidelity and converges to MLE's fidelity for 
a higher number of qubits. 



n 


MLE 


Iteration Time 


Q & D 


FP 


2 


40± 10 


1.0 ±0.5 


0.19 ±0.05 


0.17 ±0.04 


3 


1000 ± 300 


12 ±8 


0.29 ±0.03 


0.26 ±0.02 


4 


(16 ±3) 103 


180 ± 70 


1.1 ± 0.1 


1.1 ± 0.2 


5 


(35 ± 6) 104 


(4 ± 1) 103 


9.4 ±4 


12 ±4 


6 


(6.0 ±0.3) 106 


(60 ± 16) 103 


58 ±2.5 


60 ±5.4 



Table II: Runtimes in milliseconds for a state at r « 0.5 and 
Siincar ~ 0. "Iteration Time" measures how long the line 
search routine takes for each iteration of the BFGS routine. 
Abbreviations: MLE: Complete Maximum Likelihood recon- 
struction; QD: "Quick and Dirty" method; FP: "Forced Pu- 
rity". 
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n 


MLE 


Iteration Time 


Q & D 


FP 


2 


41 ± 10 


0.77 ±0.4 


0.060 ±0.005 


0.038 ± 0.006 


3 


1200 ± 200 


12 ±20 


0.17 ±0.003 


0.36 ±0.6 


4 


(19 ± 1) 103 


200 ± 100 


1.3 ±0.4 


1.1±0.1 


5 


(29 ± 2) 104 


2900 ± 900 


7.6 ±0.6 


8.8 ±2 


6 


(6.7 ±0.4) 106 


(7 ±2) 104 


65 ±5 


67 ±5 



Table III: Runtimes in milliseconds for a state at r ~ 0.5 and 
"Siincar to that along the Werner state line; abbreviations same 
as Table II. 
















1 


2 qubits 






- - - 3 qubits 






4 qubits 






- - 5 qubits 






6 qubits 






7 qubits 



0.2 0.4 0.6 0.8 




Number of Qubits 

Figure 7: For each number of qubits, this plot shows the 
number of times each projection measurement has to be re- 
peated for a state, in order to obtain an accurate estimate 
of the measurement outcome. With a 5% state error, this 
number of measurements will allow the Forced Purity routine 
to estimate the state with 90% fidelity. There is an expo- 
nential increase in how many times the experiment has to be 
repeated. 



Figure 9: Using the same procedure as in Fig. [8] we then 
performed "Forced Purity" on each state. We said that for 
the Werner state line "Forced Purity" routine improves overall 
as the number of qubits increases. This is not the case for 
pure states, and as we can clearly see "Forced Purity" slightly 
decoheres, but overall still remains at over 90% fidelity even 
for 7 qubits. 



VI. CONCLUSION 

We have demonstrated that if the experiments can 
be performed a sufficient number of times, then using 
"Forced Purity" routine, tomography can be performed 
in a quick and robust manner. However, as the entropy 
of a state increases, a much more expensive MLE rou- 
tine has to be used to perform tomography, which does 
not scale well as the number of qubits increases. Quan- 
tum computing requires only pure state tomography, for 
which we have obtained a scalable and efficient routine 
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Number of Qubits 

Figure 8: This graph shows the average fidelities of estimated 
pure states recovered using the "Quick and Dirty" approach 
for different numbers of qubits at 5% state error. In total 
11 tangle values equally spaced between and 1 were used, 
and 10 recoveries performed for each tangle value. This shows 
that while "Quick and Dirty" appears to work for 2 qubits, 
in the long run it linearly worsens and cannot be used as a 
suitable tomography algorithm. 



Appendix A: MATRIX DIFFERENTIAL 
CALCULUS THEOREMS 



The following theorems were used to derive an analytic 
expression to the gradient of the MLE function, which is 
more computationally efficient than the finite-difference 
gradient computation as the MLE function itself is 
expensive to evaluate. 
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If M is a real- valued matrix and T 
then [mini [13]: 



dM 



dM 
~dX 



+ i 



dM 

ay ■ 



X + i-Y,TeC, 



(Al) 



The following are defined for real square matrices [13j : 
dTr{X'^Y} 



dX 
dTrjX^X} 

ax 

dTriKX^Y} 
dX 

dTr{KX^X} 
dX 



= 

= 2X, 

= YK, 

= XK + XK^ 



(A2) 
(A3) 
(A4) 
(A5) 



4"-l 



We begin by observing that 

Tr{T\£)Tm= E*- 

i/=0 

which immediately implies that 
9Tr{Tt(i)r(i)} 



2U 



Hence, using matrix calculus, we have the result 
dTr{T^i)T{i)} 



dT 



2X + i- 2Y. 



(A7) 



Eq. (A7l is a compact means of stating the result of 
Eq. (A6l using a 2" x 2" matrix, where the value of the 



derivative is stored in the original position of t^, in the 
Cholesky-decomposed matrix T{t). This is the general 
idea behind all matrix calculus results we have used. We 
could have also obtained the same result by applying ma- 
trix calculus directly. For example, denote 

$ = T\t)T{i). 



Then using Eq. (|A2| and Eq. ( |A3 I 
5rr{$} 



dX 



and 



dY 



2X + i-Y -i-Y = 2X 



= i- X -i- X + 2Y = 2Y. 



This is consistent with our result in Eq. ( A7 1 . 



Next, we set out to compute '^^^i^"'^^ . Recall that 
til, ~ Ki, + i ■ Ay, hence. 



n^* = [Ki, + i ■ K){X'^ X + i ■ X'^Y ~ i ■Y'^ X + Y'^Y) 
= K^X'^X + k^y'^y + ky'^x - KX'^Y + 
i ■ {K^X'^Y + Kx'^x + ky'^y - k^y'^x). 

(A8) 



(A6) 

Applying Eq. (|A4|) and Eq. \h5\ to the real part of 



Eq. ( A8 1 we obtain 



dX 



= XK, + XKl - YK + YAl 



(A9) 



and 



dTr{Ili,<i>} 
dY 



= XA, - XAl + YK„ + YKl. (AlO) 



Observe that Vi', Iljy — yields A^ = — A^ and K,j 
, therefore 



dX 

W 



2XKy - 2YK,,, (All) 
= 2XK + 2YK„. (A12) 



Substituting the above two equations into Eq. (All we 
obtain the result described in Sec. IIIIB~2l 
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